###########################
# Table 2 and Table A5
###########################

rm(list=ls())       # clear objects in memory

#sink("~/Dropbox/COVID Peru/08_replication/03_randomization_inference_health.txt")

library(ri)         # load the RI package
library(foreign)    # package allows R to read Stata datasets

set.seed(123)

d = read.dta("clean_survey_experiment_covid_peru.dta")
d_sal = d[!is.na(d$t_ind_salud),]
d_sal2 = d_sal[!is.na(d_sal$acuerdo_cuarentena),]

table(d_sal$t_ind_salud)
table(d_sal2$t_ind_salud)

#######################################
# Health treatment Binary outcome
#######################################

Z <-  d_sal$t_ind_salud
Y <- d_sal$acuerdo_cuarentena_binary

probs <- genprobexact(Z)  
table(probs)

ate <- estate(Y,Z,prob=probs)      # estimate the ATE
ate

perms <- genperms(Z,maxiter=10000)   # set the number of simulated random assignments

Ys <- genouts(Y,Z,ate=0)    # create potential outcomes under the sharp null of no effect for any unit

distout <- gendist(Ys,perms,prob=probs)  # generate the sampling distribution  based on the schedule of potential outcomes implied by the null hypothesis

dispdist(distout,ate)       # display p-values, 95% confidence interval, standard error under the null, and graph the sampling distribution under the null

round(mean(Y[Z==1]),3)
round(mean(Y[Z==0]),3)
ate

#######################################
# Health treatment Ordinal outcome
#######################################

Z <-  d_sal2$t_ind_salud
Y <- d_sal2$acuerdo_cuarentena

probs <- genprobexact(Z)  
table(probs)

ate <- estate(Y,Z,prob=probs)      # estimate the ATE
ate

perms <- genperms(Z,maxiter=10000)   # set the number of simulated random assignments

Ys <- genouts(Y,Z,ate=0)    # create potential outcomes under the sharp null of no effect for any unit

distout <- gendist(Ys,perms,prob=probs)  # generate the sampling distribution  based on the schedule of potential outcomes implied by the null hypothesis

dispdist(distout,ate)       # display p-values, 95% confidence interval, standard error under the null, and graph the sampling distribution under the null

round(mean(Y[Z==1]),3)
round(mean(Y[Z==0]),3)
ate

sink()

